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The  present  effort  has  addressed  the  issue  of  stabilization  of  hypersonic  boundary  layer  flow  on 
account  of  porous  coating  of  the  wall,  a  phenomenon  discovered  by  Fedorov  &  Malmuth  (2001)  and 
demonstrated  experimentally  by  Rasheed  et  al.  (2001).  Two  aspects  of  the  problem  have  been  considered 
numerically,  recovery  of  two-dimensional,  essentially  nonparallel  basic  states  and  three-dimensional 
BiGlobal  instability  analysis  of  such  basic  states.  A  set  of  two-dimensional  direct  numerical  simulations 
have  addressed  slow  flow  inside  a  single  or  a  pair  of  microcavities  embedded  in  the  wall.  An  efficient 
algorithm  has  been  presented  for  the  recovery  of  basic  states  in  isolated  cavities,  driven  by  constant 
shear,  with  or  without  superimposed  small-amplitude  perturbations  originating  in  the  boundary  layer. 
Further,  the  essential  modification  of  the  boundary  layer  flow  in  the  presence  of  a  row  of  microcavities 
has  been  documented  using  spectral-element  DNS.  Parameter  ranges  have  been  identified  within  which 
it  is  permissible  to  employ  simplified  models  to  describe  flow  in  a  single  microcavity  and  calculate  the 
effect  on  the  boundary  layer  flow  in  an  analytic,  integral  manner.  Finally,  BiGlobal  instability  analysis 
has  revealed  analogies  in  the  most  unstable  eigenmodes  of  the  shear-driven  and  the  lid-driven  cavity 
flow  (Theofilis,  2000a).  Research  directions  in  the  continuing  quest  to  guide  associated  experimental 
efforts  to  control  hypersonic  boundary  layer  flow  have  been  proposed. 
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V.  Theofilis 

1.  Introduction/Theoretical  aspects 


The  motivation  for  the  present  study  is  in  the  development  of  a  better  understanding  of  the  flow  inside 
porous  coatings  which  have  recently  been  shown  to  have  very  beneficial  effects  in  controlling  hyper¬ 
sonic  boundary-layer  instabilities  (Fedorov  &  Malmuth,  2001;  Rasheed  et  al.,  2001).  In  the  present 
model,  two-dimensional  cavities  are  considered  and  a  non-parallel  (global)  linear  instability  analysis 
of  unsteady  three-dimensional  disturbances  developing  upon  the  laminar  steady  flow  inside  an  isolated 
cavity  is  performed.  The  analysis  is  distinct  from  the  plethora  of  approaches  which  monitor  the  insta¬ 
bility  of  the  shear-layer  formed  at  the  forward  lip  of  the  cavity,  and  is  akin  to  spatial  direct  numerical 
simulation  in  that  two  spatial  directions,  those  which  define  a  two-dimensional  cavity,  are  resolved 
while  the  third  is  taken  as  homogeneous  and  resolved  by  a  Fourier  expansion.  The  ability  to  resolve 
two  spatial  dimensions  in  the  global  eigenproblem  permits  a  proper  and  rational  treatment  of  factors 
affecting  flow  instability,  such  as  variable  aspect  ratio,  different  heights  of  forward/rearward  wall  of  the 
cavity  and  spacing  between  cavities.  This  approach  is  parallel  to  the  effort  of  Duck  (2002)  who  monitors 
the  external  flowfield  inside  the  boundary  layer;  the  link  between  the  two  approaches  is  provided  by 
the  boundary  conditions  at  the  open  end  of  the  cavities.  The  objective  of  the  analyses  is  to  determine 
global  three-dimensional  linear  perturbations  and  the  frequencies  of  such  disturbances  in  a  quest  to 
guide  associated  experimental  efforts  to  control  the  entire  flowfield. 

Interest  in  BiGlobal  instability  of  this  boundary  layer  has  been  raised  by  the  discovery  of  Fedorov  & 
Malmuth  (2001)  and  Rasheed  et  al.  (2001),  who  have  recently  demonstrated  theoretically  and  experi¬ 
mentally,  respectively,  that  mode  II  instability,  which  prevails  in  hypersonic  boundary-layer  flow  under 
a  variety  of  environmental  conditions,  can  be  effectively  controlled  passively  by  coating  the  surface  on 
which  the  instability  develops  by  a  porous  material.  The  objective  of  the  present  effort  is  to  provide 
detailed  numerical  predictions  of  the  flowfield  in  the  neighbourhood  of  this  porous  material.  The  key 
feature  of  this  technology  is  the  large  disparity  of  scales  between  a  typical  mode  II  wavelength  and  the 
diameter  of  each  of  the  microscopic  cavities  on  the  coating,  schematically  depicted  in  Figure  1.  The 
porous  material  can  thus  be  modelled  by  a  row  of  cavities  which  are  embedded  in  the  wall-coating;  the 
cavities  are  of  small  size  compared  with  the  thickness  of  the  boundary  layer  and  are  referred  to  as  mi¬ 
crocavities.  The  optimal  distribution  of  the  microcavities  on  the  coating  surface  is  currently  limited  by 
the  unknown  nature  of  the  interaction  between  the  flowfields  in  the  immediate  vicinity  of  neighbouring 
microcavities.  The  present  effort  addresses  this  issue  by  monitoring  separately  the  flow  regimes  firstly 
in  the  boundary  layer  over  the  porous  wall  and  secondly  inside  the  microcavity  with  the  aim  to  arrive 
at  theoretically  founded  predictions  on  control  of  mode-II-driven  laminar-turbulent  transition  in  this 
configuration.  The  boundary  condition  at  the  wall  of  the  boundary  layer  is  used  to  provide  the  link 
between  the  two  flow  regimes.  The  flow  inside  the  boundary  layer  and  its  instability  characteristics  are 
discussed  in  the  related  effort  of  Duck  (2002);  he  argues  firstly  that,  the  flow  inside  the  microcavity  is 
expected  to  be  of  a  viscous  nature,  driven  by  a  constant  shear  at  the  open  end  of  the  microcavity  and, 
secondly,  that  an  incompressible  model  may  be  used  for  the  description  of  the  basic  state  inside  the 
microcavities.  The  focus  of  the  present  numerical  effort  is  on  the  description  of  the  flowfield  inside  a 
single  microcavity  in  isolation  as  well  as  in  the  basic  flow  in  the  near-wall  region  of  a  boundary  layer 
in  which  the  minimum  nontrivial  number  of  two  microcavities  is  embedded  into  the  wall. 

Specifically,  a  number  of  simplifications  can  be  employed  in  order  to  gain  first  insights  into  the 
configuration  at  hand.  Those  pertaining  to  the  external  boundary  layer  flow  are  discussed  by  Duck 
(2002);  here  we  concentrate  on  the  simplifications  of  the  flowfield  inside  the  microcavities.  These  amount 
to  neglecting  the  three-dimensionality  of  the  geometry  for  the  purposes  of  the  present  effort  as  far  as 
the  basic  state  is  concerned.  Regarding  global  instabilities,  a  Fourier  decomposition  of  the  spanwise 
spatial  direction  is  considered.  As  has  been  mentioned,  the  flow  inside  the  microcavity  is  essentially 
incompressible  driven  by  constant  shear,  which  leaves  two  parameters  to  be  considered  besides  Reynolds 
number,  namely  the  length  over  the  depth  of  each  microcavity  and  the  spacing  between  two  neighbouring 
microcavities.  In  the  absence  of  prior  information  on  the  flow  at  hand,  here  we  focus  on  the  basic  state 
inside  a  single  microcavity  and  its  global  instability.  The  latter  is  contrasted  with  that  in  the  classic 
lid-driven  cavity  (Theofilis,  2000  a). 

The  basic  flow  is  calculated  by  numerical  solution  of  the  vorticity  transport  equation  alongside  the 
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relation  between  streamfunction  and  vorticity, 

^  r^dC_dV><9C'l  _ 

dt  \dy  dx  dx  dy  J  ^ 

=  0, 


(1.1) 

(1.2) 


where  K,  =  (l/i?e)  (T>2  +  V2) ,  while  (  =  — du/dy  +  dv/dx  is  the  vorticity  of  the  basic  flow  and  ip  is  its 
streamfunction,  for  which  u  =  dip/dy  and  v  =  — dip/dx  holds.  An  efficient  direct  numerical  simulation 
algorithm  for  the  solution  of  (1.1-1. 2)  is  presented  in  the  next  chapter. 

A  constant  shear  is  imposed  at  the  open  lip  of  a  single  microcavity  and  two  distinct  situations  are 
considered,  one  in  which  the  equations  of  motion  (1.1-1. 2)  are  integrated  in  time  until  convergence  to 
a  steady  state,  and  one  in  which  such  a  state  is  perturbed  harmonically  on  account  of  linear  instability 
originating  in  the  boundary  layer.  However,  both  basic  flows  obtained  in  this  manner  suffer  from  the 
restrictive  assumption  of  prescribed  (and  limited)  communication  between  flow  in  the  boundary  layer 
and  the  microcavity.  Furthermore,  interaction  between  the  microcavities,  which  represents  one  of  the 
main  interests  of  the  present  research,  cannot  be  addressed  by  monitoring  microcavities  in  isolation. 
The  influence  of  variation  of  geometric  parameters  of  individual  microcavities  on  the  induced  flowfield 
inside  the  boundary  layer  (and,  consequently,  on  the  stability  characteristics  of  the  flow)  is  a  further 
issue  worthy  of  investigation.  Finally,  despite  the  low  Reynolds  numbers  encountered,  it  is  by  no  means 
clear  that  the  flow  in  the  neighbourhood  of  the  microcavities  remains  steady  in  the  parameter  range  of 
interest.  In  order  to  address  all  these  issues,  basic  states  have  also  been  obtained,  considering  the  entire 
domain  encompassing  the  near-wall  portion  of  the  boundary  layer  and  a  nontrivial  minimum  of  two 
microcavities  embedded  in  the  wall. 

BiGlobal  instability  analysis  of  the  flow  is  addressed  by  numerical  solution  of  the  linear  eigenvalue 
problem  (Theofilis,  2002)  in  which  the  linear  disturbance  equations  can  be  recast  after  substitution  of 
a  decomposition  into  the  two-dimensional  0(1)  basic  quantities  q  =  (u,  v,0,p)T  and  three-dimensional 
0(e)  disturbance  quantities  q  =  q(a;,j/)e1  with  q  =  ( u,v,w,p)T  into  the  equations  of  motion, 

subtraction  of  the  basic  flow  terms  and  neglecting  terms  at  0(e2).  In  the  present  temporal  framework 
/ 3  is  taken  to  be  a  real  wavenumber  parameter  describing  an  eigenmode  in  the  2— direction  and  solve  for 
the  complex  eigenvalue  Cl.  ClT  =  5ft{0},  the  real  part  of  the  eigenvalue  is  related  with  the  frequency  of  the 
global  eigenmode  while  the  imaginary  part  is  its  growth/damping  rate;  a  positive  value  of  Cl,  =  3{f2} 
indicates  exponential  growth  of  the  instability  mode  q  in  time  t  while  Cl,  <  0  denotes  decay  of  q  in 
time.  When  the  wavenumber  vector  (3ez  is  perpendicular  to  the  plane  on  which  the  basic  flow  (u,  v,  0,p) 
develops,  as  the  case  is  in  the  present  basic  flow,  it  is  possible  to  simplify  the  two-dimensional  eigenvalue 
problem  by  rewriting  it  as  one  with  real  coefficients,  which  requires  half  the  amount  of  storage  compared 
with  the  standard  form  (Theofilis  et  al.,  2001).  The  absence  of  a  basic  flow  2;— velocity  component  in 
the  linear  operator  in  conjunction  of  the  redefinitions 


Cl  =  i  Cl,  (1.3) 

w  =  i  w  (1-4) 


result  in  the  following  generalised  real  nonsymmetric  partial  derivative  EVP  after  linearisation  and 


subtraction  of  the  basic-flow  related  terms. 

[A4  —  ( Vxu)\  u  —  ( Vyu)v  —  Uxp  =  Cl  u,  (1-5) 

—  (Vxv)u  +  [A4  —  (T>yv)\  v  —  Vyp  =  Civ,  (1-6) 

+Mw  —  (3p  =  Cl  w,  (1.7) 

T>xu  +  Vyv  —  /3w  =  0,  (1-8) 

where 

M  =  (1/Re)  (V2X  +  V2y  +  (32)  -  uDx  -  vVy.  (1.9) 


The  ability  to  solve  a  real  eigenvalue  problem  is  essential  in  the  case  of  the  present  high  Reynolds 
numbers  flow.  This  point  has  been  clearly  manifested  in  the  literature  in  the  difficulties  encountered  by 
the  investigators  who  tackled  the  problem  of  linear  instability  in  the  classic  lid-driven  cavity  to  produce 
consistent  results. 
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2.  Numerical  methods 

In  a  manner  analogous  to  classic  linear  theory,  in  which  one  spatial  direction  is  resolved  and  two  are 
considered  homogeneous  and  represented  by  Fourier  expansions,  the  details  of  the  basic  state  of  a  global 
linear  instability  analysis  crucially  condition  the  quality  of  the  stability  analysis  results.  For  this  reason 
effort  has  been  paid  to  develop  a  novel  accurate  and  efficient  numerical  approach  for  the  calculation  of 
the  basic  state,  details  of  which  are  presented  here. 

The  principles  of  the  algorithm  can  be  applied  to  recover  three  two-dimensional  velocity  compo¬ 
nents  or  indeed  a  three-dimensional  flowfield  (Ku  et  al.,  1987);  however,  for  simplicity  we  confine  the 
present  discussion  to  solutions  of  the  system  (1.1-1. 2)  which  deliver  the  basic  flowfield,  q,  inside  a  two- 
dimensional  rectangular  microcavity.  The  main  advantage  of  the  velocity-vorticity  formulation  is  that 
the  continuity  equation  is  exactly  satisfied.  However,  the  problem  of  imposition  of  boundary  conditions 
within  the  framework  of  an  overall  efficient  numerical  solution  algorithm  remains.  This  is  compounded 
by  the  fact  that  the  number  of  points  discretising  the  two  spatial  directions  in  the  subsequent  analysis 
cannot  be  increased  at  will;  while  interpolation  of  a  basic  flow  solution  obtained  on  a  very  large  number 
of  points  onto  a  modest  EVP  grid  is  one  possible  option,  it  is  more  elegant  to  avoid  the  interpolation 
procedure  altogether  and  seek  an  accurate  basic  flow  solution  on  the  same  small  number  of  discretisation 
points  on  which  the  subsequent  global  instability  analysis  is  to  be  performed.  This  appears  tailor-made 
for  a  spectral  numerical  solution  approach  (Canuto  et  al .,  1987);  in  what  follows  we  discuss  a  different 
solution  based  on  an  efficient  real-space  eigenvalue  decomposition  (EVD)  algorithm. 

Chebyshev  polynomials  have  almost  exclusively  been  used  in  the  past  in  the  context  of  spectral 
simulations  of  the  time-accurate  Navier-Stokes  and  continuity  equations,  mainly  due  to  the  availability 
of  fast  transform  algorithms  necessary  for  efficient  time-integration.  However,  for  the  present  problems 
we  have  not  restricted  ourselves  to  this  class  of  orthogonal  polynomials.  Considerable  freedom  exists  in 
the  choice  of  the  expansion  functions  and  the  associated  collocation  grids  by  using  Jacobi  polynomials 
p(i’r>  for  the  discretization  of  both  spatial  directions;  of  course,  q  =  r  =  —0.5  may  be  related  to  the 
Chebyshev-  while  q  =  r  =  0  are  the  Legendre  polynomials.  Collocation  derivative  matrices  for  both 
Jacobi  Gauss-Lobatto  and  equidistant  grids  can  be  constructed  from  first  principles;  if  (xj,  j  =  0,  •  •  •  ,  n) 
is  the  collocation  grid  chosen,  the  entries  dij  of  the  (0  :  n)  x  (0  :  n)  first-order  derivative  collocation 
matrix  T> ,  derived  analytically  from  the  interpolating  polynomial,  are 

n 

n  (xi  -  xk) 

dij  =  - — — - ,  i,j,k  =  0,---  ,n,  i^j^k,  (2.1) 

{xi  -  Xj )  II  {xj  -  xk) 

k= 0 

da  —  —  ,  i,  k  0 ,  *  *  *  ,  ^,  i  ^  L  (2.2) 

E  (Xi  -  xk) 

k= 0 


These  formulae  result  in  the  well-known  ones  if  the  analytically  known  Chebyshev  Gauss-Lobatto  grid 
(xj  =  cos  jn/n,  j  =  0 ,  ,n)  is  used  (Canuto  et  al.,  1987).  Values  of  order  m  derivatives  on  the 

collocation  grid  Xj  are  obtained  by  ( V)m . 

As  far  as  temporal  discretization  of  (1.1)  is  concerned,  the  viscous  nature  of  the  problems  in  which 
we  are  interested  introduces  scales  which  dictate  an  implicit  treatment  of  the  linear  term  in  this  equa¬ 
tion;  the  nonlinear  term  may  be  treated  explicitly.  Within  the  framework  of  solution  methods  which  do 
not  resort  to  splitting  and  introduction  of  intermediate  fields  but  rather  address  the  governing  equa¬ 
tions  directly  the  combination  of  Crank-Nicholson  (CN)  with  second-order  Adams-Bashforth  (AB2) 
or  Runge-Kutta  (RK)  schemes  has  been  extensively  used  for  the  time-integration  of  the  viscous  and 
the  convective  terms,  respectively  (Canuto  et  al,  1987).  However,  the  family  of  compact  schemes  pro¬ 
posed  by  P.  R.  Spalart  et  al.  (1991)  (SMR)  presents  more  accurate  and  more  stable  alternatives  to  the 
CN-AB2  algorithm  although  it  does  not  require  additional  computational  effort  to  the  latter  scheme. 
The  SMR  algorithm  may  be  written  in  compact  form  as 

q  =q  +  A t^Cffiq  +  A q  )  +  [iAf(q  )  +  uAf(q  )  j ,  (2.3) 
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where  the  superscript  denotes  fractional  time-step,  C{q)  and  A f(q)  are,  respectively,  the  linear  and 
nonlinear  operators  in  the  problem  to  be  solved  and  At  is  the  time-step.  The  rationale  behind  the 
derivation  as  well  as  sample  values  of  the  constants  k,  A,  p  and  v  of  a  self-starting  SMR  algorithm  may 
be  found  in  P.  R.  Spalart  et  al.  (1991).  The  time-discretization  of  (1.1)  using  (2.3)  delivers  the  following 
problem  to  be  solved  for  (£,  if)  at  each  fractional  time-step 


M  i  C'"  =  R ,  (2.4) 

M2  if'"  +  C"'  =  0,  (2.5) 

where  M\  =  d2 /dx2  +  d2/dy2  —  Re/(XAt)  and  M.2  =  d2/dx 2  +  d2/dy2,  subject  to  the  boundary 
conditions  appropriate  to  the  problem  in  consideration.  R  comprises  the  nonlinear  and  the  terms  arising 
from  the  discretization  at  previous  fractional  time-steps, 


R 


-(k/X) 


‘ d 
dx2 


dy 2 


Re/(nAt) 


c" 


(pRe/X)  (fjy(x  -  tfx'Cy)  +  ( vRejX )  ( \f'yCx 


(2.6) 


The  accuracy  of  the  overall  procedure  clearly  depends  on  the  scheme  utilised  for  calculation  of  the 
spatial  derivatives.  The  spectral  discretisation  chosen  introduces  dense  matrices  and  can  only  become 
competitive  against  other  numerical  approaches  from  the  point  of  view  of  efficiency  on  account  of 
the  existence  of  a  fast  algorithm  for  the  inversion  of  the  implicit  operators  A4i  in  (2.4)  and  A42  in 
(2.5).  While  A42  is  time-independent,  the  first  implicit  operator  is  a  function  of  a  time-step  At  that  is 
controlled  by  the  Courant-Friedrich-Lewis  condition  and  needs  to  be  inverted  at  every  time-step.  If  one 
sacrifices  the  advantage  of  an  adjustable  time-step  and  keeps  At  fixed  at  a  slightly  lower  than  its  optimal 
value,  a  powerful  EVD  algorithm  may  be  constructed  for  the  efficient  solution  of  the  incompressible 
two-dimensional  Navier-Stokes  and  continuity  equations  in  streamfunction-vorticity  formulation. 

Key  papers  on  EVD  algorithms  are  the  work  of  Haidvogel  &  Zang  (1979)  and  that  of  Ku  et  al. 
(1987).  The  first  authors  discuss  the  solution  of  Poisson’s  equation  subject  to  homogeneous  Dirichlet 
boundary  conditions  in  transform  space  while  the  second  authors  present  an  eigenvalue  decomposition 
algorithm  for  the  Poisson  equation  resulting  from  a  time— splitting  of  the  incompressible  Navier-Stokes 
and  continuity  equations  in  primitive-variables  formulation  in  the  context  of  real-space  spectral  collo¬ 
cation  using  Neumann  boundary  data.  Here  we  present  a  variant  of  the  EVD  algorithm  for  the  solution 
of  the  equations  of  motion  in  real-space  using  the  streamfunction-vorticity  formulation.  Using  this  for¬ 
mulation  within  a  direct,  as  opposed  to  time-splitting  or  iterative,  time-integration  methodology  one 
Poisson  and  one  Helmholtz  problem  are  to  be  solved  within  each  fractional  time-step.  This  minimizes 
the  necessary  computing  effort  and  makes  the  present  algorithm  one  viable  candidate  to  obtain  the 
desired  basic  states. 


(a)  Single  cavity,  steady  runs 

Physical  boundary  conditions  can  only  be  provided  for  the  stream-function  if  itself  and  its  derivatives 
on  the  domain  boundary.  We  refrain  from  discussion  of  the  Ku  et  al.  (1987)  algorithm  for  the  standard 
testbed  lid-driven  cavity  problem,  which  we  term  EVD2  for  reasons  which  will  become  apparent  in 
what  follows,  and  concentrate  on  our  novel  EVD4  algorithm  which  permits  addressing  problems  in 
which  both  Dirichlet  and  Neumann  boundary  conditions  are  imposed  on  the  stream-function  while  no 
boundary  data  are  necessary  for  the  flow  vorticity.  We  take  a  rectangular  cavity  to  be  defined  in  the 
two-dimensional  domain  ( Xi  €  [0,HR],  i  =  0,  •  •  •  ,  m)  x  (yj  G  [0,1],  j  =  0,  •  •  •  ,  n),  where  AR  is  the 
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aspect  ratio  of  the  cavity.  The  boundary  conditions  on  ip  are 


—  IpiO  —  Ipmj 

=  0, 

(2.7) 

{d2ip/dy2)in 

=  F(xi), 

(2.8) 

( dip/dx)0J 

=  0, 

(2.9) 

{dip/dy)i0 

=  0, 

(2.10) 

(dip/dx)mj 

=  0, 

(2.11) 

where  ftj  =  f(xi,yj)  represents  either  of  ip  or  £  grid-values  at  a  fractional  time-step  and  F(x)  is  a 
function  used  to  distinguish  between  a  constant-shear  boundary  condition,  F(x)  =  1,  and  a  function 
derived  from  the  external  boundary  layer  flow.  The  numerical  discretisation  of  (2. 4-2. 5)  leads  to  a  system 
of  simultaneous  equations  of  the  type 


Mf  +  /N  +  cl/  =  g,  (2.12) 

where  M  represents  the  (0  :  m )  x  (0  :  m)  discrete  analogon  of  D2,  N  represents  the  transpose  of  the 
(0  :  n)  x  (0  :  n)  discrete  analogon  of  V2,  I  is  the  identity  matrix,  c  =  —Re/(XAt)  and  g  =  R  if  (2.4) 
is  being  solved,  while  c  =  0,  g  =  — £  in  the  case  of  (2.5).  On  account  of  the  homogeneous  boundary 
conditions  (2.7)  on  ip  (2.5)  becomes 

m— 2  n— 2 

^  ^  Mki'lfiil  ^  ^  Nijtykj  —  C kl  Mkl^U  Mfom—llpm— 11  ■^ln—l'^kn—1-  (2.13) 

i= 2  j= 2 

The  boundary  conditions  (2.8-2.11)  may  be  expressed  using  the  discrete  analoga  X,Y  and  Y2  of  the 
collocation  derivative  matrices  T>x,T>y,  and  T>2,  respectively,  as  given  by  (2. 1-2. 2).  It  follows  that 


Ipll  =  ^  '  ^li^Pil  5  1pm  — 11  —  ^  '  Pm—lilpil') 

i—2  i—2 

n—2  n—2 

'iftkl  =  “1“  ^  ^  ^1  j^Pkj  i  'Ipkn—  1  =  ^ k  "t~  ^  ^  ^n—  lj'&kji 

3  = 2  3= 2 

where  Rk  and  \k  are  used  to  impose  the  boundary  condition  on  the  lid, 

—  T’fc^On-l  T  FkY0 1 


Kk  = 


A  k  = 


V  V2  V2  V  ’  K  V  V2  V2  V  ’ 

*01Inn-l  rnlron-l  *01 rnn-l  rnl10n-l 


(2.14) 

(2.15) 


(2.16) 


and  the  vectors  Jij,  Sm-u  and  ey,  en-ij  are  known  functions  of  the  entries  of  X ,  Y  and  Y2, 


i  — 


Xn 


1  X.  ■  ,  i  j  A  n 


il, 


0  i 


X01X 

mm—  1  XrnlXQrn—  1 

Aon-1^- 

= 


ToiFn2n_l  -  ^1^-1 


^m—  1  i 


Cn-lj  — 


AmlAo,  —  -Am  X 


01- 


Ami  mm—  1  -^ml^Om-1 


YAYoj  ~  YoiY^ 

Yoi ^Vi  -  ^2i^0n-i ' 


(2.17) 

(2.18) 


The  essence  of  the  EVD4  algorithm  is  to  diagonalise  the  (0  :  m  —  4) 2  matrix  M  and  the  (0  :  n  —  4) 2 
matrix  N  whose  entries  are 


Mki  =  Mki  +  MkiSu  +  Mkm-iSm-ii,  k,i  =  0,  ■■■  ,m-4,  (2.19) 

Nij  =  iVy  +  Nueij  +  Nin_itn_ij,  l,  j  =  0,  •  •  •  ,  m  —  4.  (2.20) 

The  Poisson  problem  (2.13)  becomes 

Mfki  +  fkiN  =  — C ki  -  Nn Rk  —  Nin-i Afc  =  (2-21) 
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in  which  the  non-singular  matrices  M  and  N  may  be  diagonalised  using  their  eigenvalue  decomposition 

M  =  and  N  =  {N*)v* {N*)-1 ,  (2.22) 

to  yield 

p*  (M*)-1  fki(N*)  +  (M*)~1fki(N*)  =  (M*)~1rkl(N*).  (2.23) 

As  a  consequence,  instead  of  having  to  solve  the  (to  —  3)  x  (n  —  3)  system  of  simultaneous  equations 
(2.21)  one  solves  the  (to  —  3)  x  (n  —  3)  algebraic  equations 

f*  =  r*/(p*  +  v*)  (2.24) 

for  f*  =  (M*)-1  fki{N*) ,  given  r*  =  ( M*)~1rtci(N *).  The  structure  of  the  EVD4  algorithm  is  sum¬ 
marised  in  table  2.  Clearly,  the  cost  of  this  algorithm  is  a  negligibly  small  fraction  of  the  cost  of  a  direct 
algorithm  for  the  solution  of  the  Poisson  problem.  This  is  documented  in  table  3  where  memory  and 
runtime  requirements  are  shown  for  solution  of  (1.1-1. 2)  in  the  lid-driven  cavity  problem. 

(b)  Single  cavity,  unsteady  runs 

The  algorithm  of  the  previous  section  accounts  for  flow  inside  the  cavity  set  up  by  a  constant  shear- 
driven  flow  in  the  (external)  boundary-layer.  However,  on  account  of  the  first  boundary  condition  in 
(2.7),  exchange  of  fluid  between  the  external  flow  and  the  cavity  is  prohibited.  This  can  easily  be  (partly) 
remedied  by  specifying  an  alternative  boundary  condition  at  the  North  (N— )boundary  of  the  cavity, 
which  models  unsteady  motion  generated  by  instability  in  the  external  flowfield.  One  such  model  of 
a  boundary  condition  takes  into  account  the  fact  that  linear  instability  in  the  free-stream  generates 
harmonic  motion  at  this  boundary.  Compared  with  the  previous  section  an  analogous  in  spirit,  if  not 
in  detail,  algorithm  can  be  constructed,  in  which  small-amplitude  perturbations  in  the  boundary  layer 
generate  harmonic  unsteadiness  at  the  lip  of  the  cavity.  In  turn,  such  a  motion  is  modelled  by  imposition 
of  the  inhomogeneous  boundary  conditions 


ipin  =  £COs(2nxi)cos(2nft),  (2.25) 

(d2ip/dy2)in  =  1  +  e  cos^^  +  <j>x))  cos(27t f(t  +  (ft))  (2.26) 

to  replace  the  first  boundary  condition  in  (2.7)  and  (2.8),  respectively.  Here  e  <C  1  is  a  small- 
amplitude  parameter,  /  is  an  imposed  real  frequency  of  the  motion  at  the  N—  boundary,  determined 
by  matching  with  instability  in  the  external  flow  (Duck,  2002)  and  <fx ,  (ft  are  constant  phase-shift 
parameters  in  space  and  time,  respectively. 

The  EVD4  algorithm  is  consequently  modified  in  order  to  account  for  the  boundary  conditions  (2.25- 
2.26)  while  the  last  three  of  (2.7)  and  (2.9-2.11)  are  retained.  The  matrices  M  and  N  are  then  redefined 
to  incorporate  the  unsteady  boundary  conditions,  while  M  and  N  are  identical  with  those  in  §2(a)  in 
the  limit  £  =  0.  In  this  manner  the  present  algorithm  generalises  that  of  the  previous  paragraph  without 
loss  in  efficiency  or  accuracy. 


(c)  Row  of  microcavities,  unsteady  runs 

Finally,  an  integral  part  of  the  current  research  focusses  on  a  coupled  solution  of  the  external 
boundary  layer  and  the  internal  cavity  flows.  In  this  case  a  single-domain  spectral  spatial  resolution, 
as  used  in  the  previous  sections,  is  no  longer  appropriate.  Alternative  approaches  for  the  numerical 
solution  of  such  a  problem  include  the  spectral  multidomain  algorithm  of  Theofilis  (2000a)  or  the 
spectral  element  methodologies  of  Dubois-Pelerin  (1998)  and  Karniadakis  &  Sherwin  (1999);  the  first 
two  algorithms  have  been  employed  here. 

The  minimum  nontrivial  configuration  of  two  microcavities  embedded  in  the  wall  has  been  consid¬ 
ered,  as  schematically  depicted  in  figure  7.  The  horizontal  wall  is  defined  by  y  =  0  while  microcavities 
having  length  d  and  depth  D,  placed  at  a  spacing  s  apart  have  been  considered.  Although  this  does 
not  represent  an  essential  constraint  of  the  algorithms  used,  microcavities  of  identical  dimensions  have 
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been  taken  in  order  to  reduce  the  number  of  free  parameters  in  the  problem;  in  the  same  spirit,  un¬ 
less  otherwise  stated,  square  microcavities  ( d  =  D)  have  been  considered  throughout.  In  this  manner, 
the  remaining  geometric  parameters  are  the  distance  between  the  inflow  boundary  and  the  upstream 
corner  of  the  first  microcavity,  that  between  the  downstream  corner  of  the  second  microcavity  and  the 
outflow  boundary  and,  finally,  y  =  H,  representing  the  distance  from  the  horizontal  wall  along  the 
y— direction  at  which  the  integration  domain  is  truncated.  It  is  of  interest  to  focus  future  research  in 
studying  the  modifications  to  the  solutions  obtained  herein  when  the  present  choice  of  parameters  is 
altered,  especially  by  considering  identical  cavities  of  different  dimensions  than  those  chosen  here,  or  by 
considering  microcavities  that  are  mildly  different  to  one-another,  as  the  case  is  likely  to  be  in  practical 
applications. 

Solutions  obtained  here  have  been  confined  to  two  spatial  dimensions,  but  no  assumption  of  pe¬ 
riodicity  at  the  inflow/outflow  boundaries  has  been  made;  performance  of  a  so-called  ’spatial’  direct 
numerical  simulation  (DNS)  has  been  an  explicit  requirement  in  the  early  stages  of  designing  the  current 
project  (Fedorov,  personal  communication).  Both  the  spectral  multidomain  and  the  spectral  element 
algorithm  consider  a  partition  of  the  configuration  of  figure  7  into  a  total  of  seven  parent  subdomains. 
The  first  algorithm  proceeds  by  defining  independent  conforming  two-dimensional  canonical  (global) 
mapped  Chebyshev  collocation  grids  to  discretize  each  subdomain  and  an  iterative  procedure  to  ensure 
C 1  continuity  across  domain  interfaces.  The  second  algorithm,  still  considers  the  same  space  tessellation 
but  in  addition  permits  independent  prescription  on  the  one  hand  of  the  number  of  regular  subdomains 
into  which  each  of  the  parent  subdomains  is  partitioned  and  on  the  other  hand  the  degree  of  the  polyno¬ 
mial  approximation  within  each  subdomain.  In  this  way,  the  spectral  element  approach  has  the  potential 
to  resolve  large  gradient  regions,  such  as  that  near  the  corners  at  the  lips  of  the  microcavities,  more 
efficiently  than  the  global  spectral  multidomain  approach.  In  both  cases  a  semi-explicit  time-integration 
scheme  has  been  used,  employing  low-storage  Runge-Kutta  or  backward-differences  for  the  convective 
and  Crank-Nicolson  for  the  viscous  terms.  The  boundary  conditions  used  are 


At  the  walls  : 

u  — 

v  =  0, 

(2.27) 

At  inflow  : 

u  = 

— [/ 

H  °° 

(2.28) 

V  = 

o, 

(2.29) 

At  y  =  H  : 

u  — 

Uoo, 

(2.30) 

V  = 

0, 

(2.31) 

At  outflow  : 

stress  —  free  conditions 

(2.32) 

Note  that  (2.28),  after  adjustment  of  constants,  is  consistent  with  the  constant-shear  flow  assumption 
used  throughout  this  work. 
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3.  Results 

(a)  Basic  flow 

(i)  Single  cavity,  steady  runs 

The  EVD  algorithm  discussed  in  §  2  (a)  has  been  employed,  replacing  the  uniform  tangential  velocity 
of  the  cavity  roof  (Theofilis,  2000a;  Theofilis  et  al.,  2001)  by  that  of  uniform  wall  shear  stress,  as  detailed. 
The  issue  of  the  singularity  of  the  boundary  conditions  at  the  NE  and  NW  corners  of  the  cavity  is  thus 
absent  in  the  present  calculations  and  the  time-accurate  spectral  collocation  scheme  utilized  for  the 
calculation  of  the  basic  flow  demonstrates  exponential  convergence.  Unlike  the  standard  lid-driven 
cavity,  in  which  the  condition  ipy  =  1  is  imposed  at  the  N  boundary  and  defines  the  effective  flow 
Reynolds  number,  there  is  some  ambiguity  in  the  definition  of  a  Reynolds  number  in  the  present  case. 
The  input  Reynolds  number  may  be  regarded  as 

Re  =  uv(y  =  D)d2 /v,  (3.1) 

given  that  in  our  approach  the  roof  shear  stress  is  constant  (or  ipyy  =  1);  here  v  is  the  kinematic 
viscosity  of  the  fluid.  Two  additional  Reynolds  numbers  useful  for  closer  comparisons  with  the  standard 
lid-driven  cavity  problem  may  also  be  defined,  namely  an  integral  Reynolds  number 

1  fd 

Re-mt  =  —  /  u(x,  y  =  D)  dx,  (3.2) 

v  Jx= o 

as  well  as  that  based  on  the  maximum  streamwise  velocity  component  attained  at  y  =  D, 

Re-max  =  ^  ma x{u(x,y  =  D)}.  (3.3) 

While  v  is  a  fixed  parameter  in  the  basic  flow  calculation  code,  u(x,  y  =  D )  is  unknown  a  priori  and  is 
determined  from  the  converged  basic-flow  field.  Taking  d  =  D  =  1  basic  flow  results  have  been  obtained 
using  a  rectangular  grid  comprising  upwards  of  96  collocation  points  per  spatial  direction.  These  have 
been  compared  against  the  classic  lid-driven  cavity  flow,  the  Reynolds  number  in  the  latter  case  taken 
to  be  of  O(Re-mt)  in  the  former.  As  an  example  a  comparison  between  the  constant-shear  and  lid 
driven  basic  states  in  terms  of  the  streamfunction  ip,  the  velocity  components  u  =  ipy,v  =  —  ipx  and  the 
vorticity  f  is  shown  in  figure  2.  The  shear-driven  cavity  results  have  been  generated  using  Re  =  6250, 
which  yields  an  effective  Reynolds  number  of  i?e;nt  «  300.  Both  sets  of  the  first  three  quantities  and 
the  vorticity  of  the  shear-driven  flow  are  presented  by  ten  isolines  between  zero  and  the  respective 
maxima.  The  comparison  points  to  a  general  qualitative  agreement  between  the  two  model  flows.  That 
the  input  Reynolds  number  (Re)  in  the  shear-driven  case  is  substantially  higher  than  Re-mi  shows  how 
relatively  ineffective  the  constant  shear  driving  mechanism  is.  In  anticipation  of  the  parameter  range 
necessary  to  produce  globally  unstable  flows  (Theofilis,  2000a;  Theofilis  et  al. ,  2001),  constant -shear 
basic  states  were  obtained  at  substantially  larger  values  of  the  unit  Reynolds  number,  as  shown  in 
figure  3.  A  monotonic  increase  of  both  Re-mi  and  f?emax  as  Re  =  1/v  increases  can  be  seen,  which 
enhances  the  viscous  assumptions  discussed  earlier.  The  extent  to  which  the  quantitative  differences 
between  the  basic  states  modify  the  stability  characteristics  of  the  lid-driven  cavity  flow  is  examined 
in  §3  (b). 

(ii)  Single  cavity,  unsteady  runs 

Before  doing  so,  though,  basic  flow  results  have  been  obtained  using  the  DNS  algorithm  of  §  2  (6). 
The  objective  has  been  to  verify  the  ability  of  the  code  to  couple  unsteadiness  generated  by  instability  in 
the  external  field,  which  penetrates  inside  the  microcavity,  or  unsteadiness  due  to  resonance  originating 
inside  the  microcavity  which  propagates  in  the  boundary  layer.  As  will  be  discussed  in  §  3  ( b ),  this 
coupling  is  pivotal  in  order  to  combine  the  analytic/numeric  approach  of  Duck  (2002)  in  the  external 
boundary  layer  flow  with  the  present  fully  numerical  approach  required  to  describe  flow  inside  the 
microcavity. 
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t 

£ 

/ 

<  20 

0 

0 

20  <  t  <  40 

10"4 

1 

40  <  t  <  60 

10"4 

0.5 

>  60 

0 

0 

Table  1.  Parameters  for  the  single-cavity  unsteady  simulations 

To  this  end,  a  cavity  of  aspect  ratio  (depth  over  length)  of  2  has  been  considered  and  constant 
parameter  values  <px  =  eft  =  0  have  been  used.  The  magnitude  of  the  parameter  e,  on  the  other  hand, 
is  expected  to  have  a  critical  qualitative  influence  on  the  results  of  the  simulations:  taking  a  low  value 
should  result  in  linear  response  of  the  flow  to  the  external  forcing  (see  2.25-2.26),  where  the  imposed 
frequency  /  is  in  resonance  with  that  by  which  the  cavity  responds;  progressively  increasing  e  will  result 
in  departure  from  this  behaviour  towards  nonlinearity. 

Representative  results  of  simulations  are  shown,  in  the  form  of  time-dependence  of  a  flow  quantity 
at  a  given  position  in  the  flowfield  in  figure  4  (here  u  is  monitored,  although  others  exhibit  qualitatively 
analogous  behaviour)  and  in  the  form  of  the  spatial  distribution  of  ip,  u,  v  and  £  at  two  characteristic 
times,  t  =  20  and  t  =  40,  in  figures  5  and  6,  respectively;  v  =  10-2,  i.e.  Re  =  100  using  the  definition 
(3.1),  and  the  parameters  of  table  1  have  been  used  for  these  simulations.  In  figure  4  it  can  be  seen 
that,  starting  from  rest,  at  t  <  20  a  steady  state  solution  is  approached;  components  of  this  solution 
are  shown  in  figure  5.  When  non-zero  amplitude  forcing  is  applied  at  the  N  boundary  of  the  cavity,  the 
flow  responds  in  a  periodic  manner  as  would  be  expected  from  (2.25-2.26).  After  rather  short  transients 
following  the  (Heaviside)  changes  of  £  at  t  =  20  and  of  /  at  t  =  40,  the  flow  settles  to  two  different 
harmonic  motions  with  periods  T2 o_4o  =  1//  =  1  at  20  <  t  <  40  and  T40_6o  =  1//  =  2  at  40  <  t  <  60. 
Snapshots  of  the  solution  at  t  =  40  are  shown  in  figure  6.  At  t  =  60  the  forcing  is  removed  and,  by 
the  end  of  the  simulation,  the  flow  approaches  the  same  steady-state  reached  at  t  =  20~.  At  these 
parameters,  this  steady  state  corresponds  to  Relnt  «  14. 

Using  the  same  parameter  values  but  higher  resolution  has  demonstrated  convergence  of  the  results 
presented  and  points  to  the  reliability  of  the  statements  put  forward  in  this  section.  Further,  qualitatively 
analogous  results  obtained  at  parameter  values  different  to  those  presented  in  table  1  point  to  the 
generality  of  the  conclusions  reached,  namely  that  the  imposition  of  harmonic  motion  at  the  open  end 
of  the  cavity  sets  up  a  linear  response  of  the  flow  inside  the  cavity  in  resonance  with  the  imposed 
frequency,  provided  that  the  forcing  amplitude  e  is  kept  small.  Another  potentially  significant  result  for 
subsequent  modelling  is  the  finding  that  the  bulk  of  the  fluid  flow  motion,  both  unforced  and  forced,  is 
confined  within  a  depth  approximately  equal  to  the  length  of  the  microcavity,  as  shown  in  the  results 
of  figures  5  and  6. 

(iii)  Two  microcavities,  unsteady  runs 

Unsteady  simulations  have  been  performed,  considering  the  entire  flowfield  which  encompasses  the 
near-wall  part  of  the  boundary  layer  and  two  microcavities  embedded  in  the  wall.  Such  simulations 
require  orders-of-magnitude  larger  computing  effort  compared  with  the  single-domain  calculations 
and,  hence,  a  limited  amount  of  results  monitoring  variations  of  the  essential  flow  parameters  have  been 
obtained.  Nevertheless,  these  results  suffice  to  highlight  essential  differences  between  the  numerical 
solutions  obtained  in  the  previous  two  sections  and  the  present  more  elaborate  approach  and  point  to 
the  direction  of  future  research.  Specifically,  two  questions  have  been  formulated  and  investigated: 

•  First,  how  does  the  flowfield  in  the  near-microcavity  region  compare  with  solutions  obtained 
numerically  using  the  assumptions  of  the  previous  section,  or  those  used  in  the  analytic  model  of 
Fedorov  &  Malmuth  (2001)? 
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•  Second,  how  do  the  geometric  parameters  s  (spacing  between  successive  microcavities)  and  D 
(depth  of  a  microcavity)  affect  the  flowfield  of  a  microcavity  (including  the  induced  boundary- 
layer  flow)  and  the  interaction  of  this  flowfield  with  neighbouring  microcavities? 

Answers  to  these  questions  can  provide  guidance  to  theory  and  experiment  alike.  Further  mod¬ 
elling  can  benefit  from  identification  of  the  limitations,  in  terms  of  the  minimum  spacing/depth  of 
microcavities,  for  which  the  integral  conditions  used  in  the  analysis  (Fedorov  &  Malmuth,  2001)  hold, 
before  interaction  between  the  flowfields  of  the  microcavities  invalidates  the  assumption  of  isolated 
microcavities  under  which  these  integral  conditions  are  derived.  On  the  other  hand,  experiment  can 
use  information  on  the  minimum  spacing  between  microcavities  when  manufacturing  materials  whose 
stability  properties  are  reliably  predicted  by  theory.  Clearly,  a  major  unanswered  question  in  this  con¬ 
text  is  the  influence  of  three-dimensionality  on  conclusions  reached  here.  This  is  an  issue  that  is  partly 
addressed  in  §  3  (6)  and  should  be  further  expanded  in  future  work. 

Taking  L  =  d  =  D  =  1,  some  numerical  experimentation  was  necessary  in  order  to  determine  the 
remaining  geometric  parameters,  such  that  (a)  the  location  of  the  inflow  boundary  does  not  interfere 
with  the  flowfield  in  the  neighbourhood  of  the  upstream  microcavity,  on  the  one  hand  permitting 
the  boundary  layer  to  develop  on  the  upstream  wall  of  the  configuration  and  on  the  other  hand  not 
placing  excessive  resolution  requirements  on  account  of  a  long  upstream  domain,  (6)  the  location  of  the 
outflow  boundary  be  placed  well  downstream  of  the  downstream  cavity,  such  that  the  outflow  boundary 
conditions  are  physically  plausible  and,  finally,  (c)  the  location  at  which  the  domain  is  truncated  in  the 
y— spatial  direction  be  such  that  activity  taking  place  in  the  neighbourhood  of  the  microcavities  is  not 
affected  by  this  (artificial)  domain  truncation. 

Having  ensured  all  three  requirements,  the  remaining  physical  parameters  with  which  numerical 
experiments  were  performed  were  the  Reynolds  number,  defined  as 


Re 


Uqq  d 
v 


(3.4) 


and  the  distance  s  between  the  microcavities.  Additional  numerical  parameters  that  needed  to  be 
determined  were  the  time-step  At  of  the  unsteady  simulations  and  the  number  of  subdomains  in  which 
each  of  the  parent  subdomains  (pi  —  pr  in  figure  7)  was  to  be  subdivided,  as  well  as  the  degree  of 
polynomial  approximation  within  each  subdomain.  The  last  two  parameters  determine  the  convergence 
properties  of  the  algorithm  and  the  accuracy  of  the  solutions  obtained  and  are  (alongside  the  Reynolds 
number)  intimately  linked  with  At  via  the  CFL  condition.) 

Results  have  been  obtained  for  Re  £  [10, 103]  and  s  £  [1,4]  while,  in  order  to  verify  the  conclusion 
put  forward  in  the  previous  section  regarding  the  significance  of  the  microcavity  depth  parameter, 
simulations  comparing  the  results  of  D  =  1  and  2,  keeping  all  other  parameters  identical,  have  also 
been  performed.  In  all  calculations  to  be  presented  the  spectral  element  code  has  been  used,  since  it 
has  been  found  that  this  methodology  is  substantially  more  efficient  in  resolving  sharp  gradients  (such 
as  those  at  the  lips  of  the  microcavities)  by  local  refinement  of  the  grid,  as  opposed  to  the  spectral 
multidomain  code  in  which  such  geometric  singularities  impose  global  requirements  on  the  grid  and, 
consequently,  higher-than-necessary  resolution  away  from  the  singularity.  In  this  respect  each  of  the 
7  parent  subdomains  p\  —  p-j  was  further  subdivided  into  42  subdomains,  within  each  of  which  flow 
quantities  were  resolved  using  a  7th  degree  polynomial.  These  values  represent  a  compromise  between 
accuracy  and  efficiency;  resolving  a  parent  subdomain  using  22  subdomains  and  lower  degree  polynomials 
produces  only  qualitatively  correct  results,  while  further  increasing  the  number  of  subdomains  (62  x  7  = 
252  being  the  next  possibility,  compared  with  the  currently  used  42  x  7  =  112  subdomains)  would  be 
unnecessarily  expensive.  In  the  parameter  ranges  explored  the  unsteady  algorithm  has  yielded  stationary 
solutions  only;  at  convergence  successive  time-step  results  exhibit  a  relative  time  variation  of  less  than 
0.1  x  10"6. 

Of  the  results  obtained  most  relevant  for  the  problem  at  hand,  as  has  been  stressed  here  and  by 
Duck  (2002),  are  low  Reynolds  number  calculations.  The  spatial  distribution  of  u(x,y )  and  v(x,y)  at 

f  In  the  case  of  multidomain  solutions  the  global  spectral  collocation  grid  utilised  to  resolve  each  parent  subdomain  is 
the  only  means  of  improving  the  solution  accuracy. 
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Re  =  50  and  s  =  D  =  1  is  presented  in  figure  8;  the  slow  motion  of  flow  inside  the  microcavity  can 
be  appreciated  in  this  result.  Closer  examination  of  the  normal  (to  the  wall  on  which  the  boundary 
layer  develops)  velocity  component  reveals  two  interesting  aspects  of  the  flow.  On  the  one  hand,  it 
can  be  seen  that  little  communication  takes  place  between  fluid  in  the  boundary  layer  and  that  in 
the  cavity;  fluid  in  the  microcavity  is  found  to  be  in  near-solid  body  rotation,  as  expected  by  the 
smallness  of  the  Reynolds  number.  On  the  other  hand,  compensation  of  this  phenomenon  takes  place 
in  the  boundary  layer,  where  the  presence  of  the  microcavity  appears  to  exert  influence  on  the  flow 
over  a  relatively  large  part  of  the  domain  outside  the  microcavity  itself.  The  first  result  is  in  line  with 
the  assumption  of  the  previous  two  paragraphs  and  those  of  Fedorov  &  Malmuth  (2001)  and  Duck 
(2002).  The  second  result  is  expected  to  have  a  strong  influence  on  the  stability  characteristics  of  the 
boundary  layer  and  can  only  be  predicted  numerically  by  simulations,  such  as  these  presented  herein. 
Further  inspection  of  the  normal  velocity  component  result  appears  to  suggest  that  s  =  1  is  a  distance 
at  which  the  fluid  motion  induced  in  the  boundary  layer  by  one  microcavity  interacts  nonlinearly  with 
that  induced  by  its  neighbour.  In  other  words,  placement  of  microcavities  at  distances  s  <  1  apart  will 
prohibit  applicability  to  a  row  of  microcavities  (in  an  integral  manner)  of  analytical  models  derived  on 
the  basis  of  a  single  microcavity.  Next,  the  effect  of  s  on  the  solutions  obtained  at  constant  Re  and  D 
is  examined;  representative  results  are  presented  in  figure  9  at  Re  =  20,1?  =  1.  It  can  be  seen  that 
s  >  4  is  necessary  in  order  for  the  induced  flowfield  of  the  downstream  microcavity  to  become  similar 
to  that  of  the  upstream  microcavity.  Unless  three-dimensional  instability  modifies  these  basic  flows, 
this  result  suggests  that  using  such  a  spacing  between  microcavities  can  ensure  that  the  model  used  to 
describe  the  induced  flowfield  of  a  single  microcavity  can  be  applied  in  an  integral  manner  to  describe 
the  interaction  between  the  boundary  layer  and  a  porous  material.  The  Reynolds  number  effect  on 
the  conclusions  put  forward  so  far  is  next  determined  by  obtaining  solutions  in  the  range  of  interest, 
Re  £  [10,  20],  keeping  s  =  2,  D  =  1;  the  results  at  the  extreme  Reynolds  number  values  are  presented  in 
figure  10.  It  can  be  seen  that,  at  least  in  the  parameter  ranges  investigated,  in  which  no  unsteadiness 
occurs,  the  Reynolds  number  effect  is  substantially  less  pronounced  than  that  of  the  spacing  between 
microcavities.  Finally,  figure  11  essentially  confirms  the  conclusions  of  the  previous  sections,  namely 
that  the  bulk  of  the  activity  inside  the  microcavity  takes  place  within  a  depth  equal  to  the  width  of 
the  microcavity;  quantitative  inspection  of  the  solution  in  y  £  [—2,  —1]  revealed  a  rapid  decay  of  all 
flow  quantities  with  depth.  The  results  obtained  have  been  analysed  quantitatively  in  different  ways,  f.e. 
monitoring  one-dimensional  profiles  of  velocity  components  as  a  function  of  the  wall-normal  coordinate 
y.  This  task,  however,  is  physically  relevant  only  if  three-dimensional  instability  does  not  destroy  the 
two-dimensional  basic  states  obtained;  one  means  of  addressing  this  instability  is  BiGlobal  instability 
theory  (Theofilis,  2002)  to  which  we  return  now. 


( b )  Global  linear  instability 

Three-dimensionality  can  be  introduced  in  the  results  of  the  previous  section  by  receptivity,  followed 
by  either  linear  instability  or  transient  growth.  Examination  of  all  relevant  mechanisms  is  beyond  the 
scope  of  the  present  work;  focussing  on  BiGlobal  linear  instability,  one  result  that  may  be  inferred  from 
the  steadiness  of  the  basic  states  obtained  is  that  the  two-dimensional  BiGlobal  eigenmode  (/ 3  =  0) 
is  stationary  (Theofilis,  2000  b).  On  the  other  hand,  linear  amplification  of  three-dimensional  BiGlobal 
eigenmodes  (having  (3  ^  0)  can  be  examined  by  numerical  solution  of  the  partial  derivative  eigenvalue 
problem  (1.5-1. 8). 

At  present,  only  single-domain  solutions  of  the  EVP  have  been  obtained,  using  the  steady  basic  states 
calculated  in  §  3  (a)  (i)  and  monitoring  the  relative  stability  properties  of  shear-driven  and  lid-driven 
cavity  flow.  The  straightforward  viscous  boundary  conditions  on  the  disturbance  velocity  components 
have  been  imposed  on  the  three  solid  walls,  while  the  issue  of  boundary  conditions  on  the  N— boundary 
of  the  shear-driven  cavity  is  expected  to  be  the  key  linking  the  flowfield  of  an  isolated  microcavity 
with  that  prevailing  in  the  external  flow,  studied  by  Duck  (2002).  In  the  absence  of  prior  information, 
calculations  reported  here  have  been  performed  by  imposing  homogeneous  Dirichlet  or  extrapolation 
of  the  disturbance  velocity  components  from  the  cavity  interior;  appropriate  compatibility  conditions 
for  the  disturbance  pressure  have  been  used  to  close  the  system  to  be  solved  in  both  cases.  In  the 
framework  of  a  temporal  global  instability  analysis  the  only  parameters  of  the  problem  are  the  flow 
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Reynolds  number  Re  and  the  spanwise  periodicity  length  LZl  associated  with  a  wavenumber  (3  =  2tt /Lx, 
while  the  complex  temporal  eigenvalue  Q  and  the  two-dimensional  amplitude  functions  ( u,v,w,p)T  of 
the  disturbance  quantities  are  determined  by  the  numerical  solution  of  the  partial-derivative  eigenvalue 
problem. 

At  Re idc  =  item*  ~  300  the  lid-driven  cavity  flow  is  known  to  be  stable  to  both  two-  and  all 
three-dimensional  perturbations.  Comparisons  between  the  two  model  flows  are  only  possible  in  the 
case  of  homogeneous  Dirichlet  boundary  conditions  being  imposed  at  the  N  boundary  in  the  shear- 
driven  case.  The  modification  of  the  (least  damped)  leading  eigenvalues  of  the  eigenspectrum  at  (3  =  7.5 
(approximately  corresponding  to  least-damped  conditions  for  mode  T2  (Theofilis,  2000a;  Theofilis  et  al, 
2001))  is  shown  in  the  upper  part  of  figure  12.  Owing  to  the  disparity  between  Re idc  and  Re  the  shear- 
driven  flow  experiences  much  weaker  damping  compared  with  its  lid-driven  counterpart.  Noteworthy 
is  also  the  difference  between  the  two  flows  in  terms  of  the  prevalence  of  stationary  and  low-frequency 
disturbances  (in  the  neighbourhood  of  Or  =  0)  in  the  shear-driven  cavity  case  at  Re  =  6250.  At 
Re  =  25  x  104,  on  the  other  hand,  the  flow  is  marginally  unstable  at  (3  =  7.5;  the  corresponding  spectrum 
is  shown  in  the  lower  part  of  figure  12.  The  amplitude  function  of  the  spanwise  disturbance  velocity 
component  w  and  the  disturbance  pressure  p  of  the  unstable  global  eigenmode  is  shown  in  figure  13. 
The  complexity  of  the  (periodic  in  z)  flowfield,  setup  by  global  instability  inside  the  microcavity,  can 
be  seen  in  this  figure,  while  the  analogies  of  the  structure  of  the  global  eigenmode  with  that  in  the 
lid-driven  cavity  (Theofilis,  2000a;  Theofilis  et  al. ,  2001)  can  also  be  appreciated. 

The  currently  imposed  boundary  condition  of  impermeability  at  the  roof  of  the  cavity  has  been 
extended  to  that  suggested  by  Duck  (2002), 


v  —  Np  =  0,  (3.5) 

which  can  provide  the  link  between  the  external  and  internal  disturbance  flowfields.  However,  in  view 
of  the  basic  flow  results  of  §  3  (a)  (iii),  further  stability  calculations  have  not  been  performed  until  the 
issue  of  the  appropriate  basic  state  to  be  used  for  such  stability  calculations  has  been  addressed  in  a 
satisfactory  manner,  a  task  constituting  one  of  the  possible  interesting  extensions  of  the  present  effort. 
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4.  Conclusions 

This  report  may  be  regarded  as  a  first  step  in  modelling  the  flow  inside  porous  media  coatings  which 
are  known  to  be  useful  in  controlling  hypersonic  (in  particular  mode  II)  boundary-layer  flow  instabil¬ 
ities.  The  practically  incompressible  nature  of  the  flow  inside  the  microcavities  strongly  suggests  that 
temperature  perturbations  in  the  outer  flow  are  unlikely  to  be  substantially  affected  by  the  presence  of 
the  microcavities. 

The  first  set  of  results  obtained  has  assumed  the  outer  flow  to  affect  the  inner  (cavity)  flow,  but  not 
vice-versa.  The  nature  of  two-dimensional  basic  flow  inside  a  microcavity  driven  by  a  uniform  shear- 
stress  on  one  face  was  established  by  use  of  an  accurate  and  efficient  eigenvalue  decomposition  algorithm 
(EVD)  for  direct  numerical  simulation.  Details  of  the  novel  aspects  of  the  EVD  algorithm  have  been 
presented.  Some  BiGlobal  instability  results  of  this  flow  have  also  been  obtained;  they  strongly  suggest 
the  flow  to  be  quite  stable  in  likely  practical  regimes  (Reynolds  numbers),  particularly  when  compared 
with  the  lid-driven  cavity  problem  (Theofilis,  2000a;  Theofilis  et  al. ,  2001). 

The  second  set  of  results  obtained  has  permitted  small-amplitude  harmonic  perturbations  super¬ 
imposed  upon  the  solutions  obtained  in  the  first  leg  of  the  investigations.  An  isolated  microcavity  has 
been  shown  to  respond  to  external  forcing  in  a  linear  manner,  periodic  flow  being  set  up  inside  the 
cavity  at  the  imposed  external  frequency.  This  extension  is  necessary  in  order  to  permit  a  more  inter¬ 
active  regime,  especially  with  respect  to  the  perturbed  flow.  In  particular  a  more  realistic  boundary 
condition  recently  suggested  by  Duck  (2002)  can  be  used  to  relate  the  wall-normal  component  of  the 
disturbance  velocity  and  disturbance  pressure  on  the  external  face  of  the  cavity,  as  derived  from  the 
present  numerical  approach,  with  those  of  the  external  flowfield  analysed  by  Duck  (2002). 

Finally,  the  basic  flow  results  obtained  using  the  previous  two  sets  of  approximations  have  been  put 
in  perspective  by  performing  (unsteady)  DNS  of  the  entire  flowfield,  which  encompasses  both  the  near¬ 
wall  boundary  layer  and  two  microcavities  embedded  in  the  wall.  It  has  been  shown  that  the  parameter 
on  which  the  flowfield  most  critically  depends  on  is  the  spacing  between  microcavities;  beyond  a  spacing 
s  «  4  (which  scales  with  other  geometric  parameters  of  the  flow)  it  appears  to  be  permissible  to  study 
microcavities  in  isolation  from  each  other  (but  interacting  with  the  boundary  layer)  and  calculate  the 
effect  of  the  induced  flow  motion  on  the  boundary  layer  in  an  analytical  (integral)  manner.  Other 
parameters,  such  as  the  depth  of  the  cavities  and  the  Reynolds  number  have  been  found  to  have  a  lesser 
impact  on  the  flow,  at  least  in  the  parameter  range  of  interest  where  steady  states  prevail. 

Future  work  aiming  to  aid  optimization  of  the  process  of  microcavity  scale  and  distribution  should 
address  the  key  issue  left  open  in  the  current  investigations,  namely  three-dimensionality.  This  can 
be  studied  either  as  harmonic  modification  in  the  third  spatial  direction  of  the  present  (single-  and 
multiple-microcavity  configuration)  two-dimensional  basic  states  or  through  consideration  of  essentially 
inhomogeneous  three-dimensional  microcavities  of  different  geometric  characteristics  embedded  in  the 
boundary  layer. 
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1  Pre-processing  Stage 

a.  Set  up  the  matrices  M  and  N, 
calculate  their  EVD  and  store  the  results 

b.  Set  up  the  matrices  M  and  N, 
calculate  and  store  their  EVD 

2  Time-advancement 

i.  First  fractional  time-step 

a.  Given  an  initial  ip  calculate  £  =  —  \/2ip 

or  read-in  (ip  ,  C  )  data  generated  at  an  earlier  simulation 

b.  Calculate  derivatives  of  ip  and  £  and  form  R 

c .  Use  EVD2  to  solve  V2£  =  R 

d.  Use  EVD4  to  solve  X/2ip  =  — £ 

e.  Overwrite  (ip  ,  £  )  by  (ip  ,  £  )  and  (ip  ,  C  )  by  (ip  ,  £  )  respectively 

ii.  Second  fractional  time-step 

//  //  /  / 

a.  Given  ip  ,ip  and  £  calculate  their  derivatives  and  form  R 

c .  Use  EVD2  to  solve  V2£  =  R 

d.  Use  EVD4  to  solve  V2ip  =  —  ( 

e.  Overwrite  (ip  ,  £  )  by  (ip  )  and  (ip  ,  C  )  by  (ip  )  respectively 

iii.  Third  fractional  time-step 

a.-e.  Same  as  2  ii. 

iv.  Check  convergence  in  time  and  either  Go  To  2  i.  or  Exit 

Table  2.  An  algorithm  for  the  solution  of  the  two-dimensional  incompressible  Navies-Stokes  and  continuity 
equations  using  the  streamfunction-vorticity  formulation  and  eigenvalue  decomposition 


F61775-01-WE049  Final  Report 


Numerical  prediction  of  the  hypersonic  boundary -layer  over  a  row  of  microcavities 


19 


SUN  Sparc  10  NEC  SX4 

problem  EVD  direct  inversion  EVD  direct  inversion 


size 

size 

(Mb) 

time 

(sec) 

size 

(Mb) 

time 

(sec) 

size 

(Mb) 

time 

(sec) 

size 

(Mb) 

time 

(sec) 

16  x  16 

0.4 

4.4 

i.i 

3.7 

4.03 

0.03 

4.03 

0.1 

24  x  24 

0.5 

5.4 

3.6 

10.5 

5.03 

0.14 

6.03 

0.3 

32  x  32 

0.6 

6.6 

10.0 

39.7 

5.03 

0.25 

13.03 

1.1 

48  x  48 

0.8 

15.3 

46.9 

460.6 

6.03 

0.56 

48.03 

8.4 

64  x  64 

1.1 

31.8 

143.9 

5203.5 

6.03 

1.08 

140.03 

40.7 

96  x  96 

1.8 

143.3 

(*) 

(*) 

8.03 

2.48 

680.03 

417.4 

128  xl28 

2.8 

523.1 

(*) 

(*) 

8.03 

4.41 

(*) 

(*) 

Table  3.  Comparison  of  memory  and  runtime  requirements  for  a  single  solution  of  a  two-dimensional  Poisson 
equation  using  direct  inversion  and  the  EVD±  algorithm  on  one  processor  of  a  workstation  and  a  supercomputer. 
Asterisks  denote  that  the  respective  problem  does  not  fit  in  the  available  memory  on  the  workstation  or  that  it 
cannot  be  solved  within  the  existing  batch  queue  time-limit  on  the  supercomputer. 
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Figure  2.  Comparison  of  the  basic  states  for  the  shear-driven  at  Re jnt  =  300  (upper  row)  and  the  lid-driven 
cavity  flow  at  the  same  Reynolds  number  (lower  row).  Left  to  right  column:  ip,u,v,(. 
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open  cavity  flow  driven  by  constant  shear 


Figure  4.  Periodic  flow  set  up  inside  an  aspect -ratio-two  microcavity,  as  a  consequence  of  a  shear  which  is 
harmonically-dependent  on  time;  shown  is  u(x  =  0.5,  y  =  1;  t). 
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open  cavity  flow  driven  by  constant  shear 


open  cavity  flow  driven  by  constant  shear 


Figure  5.  The  steady-state  solution  approached  at  t  =  20  .  Components  shown  are  ip  (upper  left),  u  (upper 

right),  v  (lower  left),  (  (lower  right). 
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open  cavity  flow  driven  by  constant  shear 


open  cavity  flow  driven  by  constant  shear 


open  cavity  flow  driven  by  constant  shear 


open  cavity  flow  driven  by  constant  shear 


Figure  6.  Snapshots,  at  t  =  40,  of  the  periodic  flowfield  set  up  inside  the  microcavity  by  harmonic  motion  at 
the  cavity  lip.  Components  shown  are  if  (upper  left),  u  (upper  right),  v  (lower  left),  C,  (lower  right). 
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Figure  7.  Sketch  of  the  geometry  considered. 
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Figure  8.  Re  =  50  flow  over  two  unit  square  microcavities,  set  apart  by  s  =  1.  Upper:  streamwise  velocity 

component;  lower:  normal  velocity  component. 
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Figure  9.  Effect  of  spacing  s  €  [2,  4]  on  a  two-unit-square-microcavity  configuration  at  Re  =  20. 
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Figure  10.  Reynolds  number  effect  on  the  normal  velocity  component  in  a  two-unit-square-microcavity 

configuration  at  s  =  2. 
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Figure  11.  Effect  of  depth  D  =  1  and  2  on  in  a  two-microcavity  configuration  at  Re  =  20,  s  =  2. 
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LATERAL  VELOCITY:  -0.0202572.  -000941906  0.30=41906  O.C282572 


PRESSURE:  -0.00021695+  -7.231SE-05  7.2318E-Q5  0.00021695^ 


Figure  13.  Four  isosurfaces  of  the  spanwise  disturbance  velocity  component  (upper)  and  disturbance  pressure 
(lower)  generated  by  global  instability  inside  a  square  microcavity  at  (Re  =  25  x  10 4 , /3  =  7.5).  Levels  are 
equidistributed  between  the  respective  minimum  and  maximum  values. 
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Latin  Symbols 

d,  D 

2D  microcavity  width,  depth 

f 

(dimensionless)  frequency 

H 

2D  microcavity  domain  height 

Re 

Reynolds  number 

t 

time 

x,  y 

Cartesian  coordinates/  collocation  points 

U,  V,  w 

Cartesian  velocity  components 

Greek  Symbols 

P 

Spanwise  wavenumber 

At 

Timestep 

8 

Small- amplitude  parameter 

K,  L,  (I,  V 

Time-integration  constants 

<f> 

Phase-shift  parameter 

?  ? 

•  9  * 

Stream-function,  vorticity 

a 

Complex  eigenvalue  of  the  global  eigenvalue  problem 

Qr 

Real  BiGlobal  frequency 

Qi 

Real  BiGlobal  amplification/damping  rate 

Calligraphic  symbols 

D 

(analytic)  collocation  derivative 

K,  L,  M,  N 

Linear  operators 

Superscripts 

i  99  999 

9  9 

Fractional  timestep  quantity 

- 

Basic  flow  quantity 

f  ~ 

Disturbance  flow  quantity 

* 

Eigenvalue  vector 

Subscripts 

1  i  i  k  1  m  ri  1 
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